#include <QFile>
#include <QTextStream>
#include <QtOpenGL>
#include "grid.h"

Grid::Grid(const QVector3D& center, QObject *parent) :
    Object3D(center, parent)
{
    m_iso = 0.0;
}

void Grid::loadFromFile(const QString& fn)
{
    QFile f(fn);
    if(!f.open(QIODevice::ReadOnly))
        return;

    QTextStream s(f.readAll());
    f.close();

    s >> m_tx >> m_ty >> m_tz;

    m_data = new float**[m_tx];
    for(int i = 0; i < m_tx; ++i)
    {
        m_data[i] = new float*[m_ty];
        for(int j = 0; j < m_ty; ++j)
        {
            m_data[i][j] = new float[m_tz];
            for(int k = 0; k < m_tz; ++k)
                s >> m_data[i][j][k];
        }
    }
}

void Grid::loadFromRawFile(const QString& fn)
{
    QFile f(fn);
    if(!f.open(QIODevice::ReadOnly))
        return;

    char v;
    m_tx = 256;
    m_ty = 256;
    m_tz = 256;

    m_data = new float**[m_tx];
    for(int i = 0; i < m_tx; ++i)
    {
        m_data[i] = new float*[m_ty];
        for(int j = 0; j < m_ty; ++j)
        {
            m_data[i][j] = new float[m_tz];
            for(int k = 0; k < m_tz; ++k)
            {
                QByteArray b = f.read(1);
                m_data[i][j][k] = b.toFloat();
            }
        }
    }
}

void Grid::setIso(float iso)
{
    m_iso = iso;
}

void Grid::incIso(int steps)
{
    m_iso += 0.005*steps;
}

void Grid::drawGeometry(void) const
{
    QVector3D pontos[8];
    float densidades[8];
    if(m_tx == 0 || m_ty == 0 || m_tz == 0)
        return;
    float dx = 1.0/m_tx;
    float dy = 1.0/m_ty;
    float dz = 1.0/m_tz;

    //cubo externo [-1]^3 a [1]^3
    for(int i = 0; i < m_tx-1; ++i)
    {
        float x0 = dx * 2 * i - 1.0;
        float x1 = dx * 2 * (i+1) - 1.0;
        for(int j = 0; j < m_ty-1; ++j)
        {
            float y0 = dy * 2 * j - 1.0;
            float y1 = dy * 2 * (j+1) - 1.0;
            for(int k = 0; k < m_tz-1; ++k)
            {
                float z0 = dz * 2 * k - 1.0;
                float z1 = dz * 2 * (k+1) - 1.0;

                pontos[0] = QVector3D(x0, y0, z1);
                pontos[1] = QVector3D(x1, y0, z1);
                pontos[2] = QVector3D(x0, y1, z1);
                pontos[3] = QVector3D(x1, y1, z1);
                pontos[4] = QVector3D(x0, y0, z0);
                pontos[5] = QVector3D(x1, y0, z0);
                pontos[6] = QVector3D(x0, y1, z0);
                pontos[7] = QVector3D(x1, y1, z0);

                densidades[0] = m_data[i]  [j]  [k+1];
                densidades[1] = m_data[i+1][j]  [k+1];
                densidades[2] = m_data[i]  [j+1][k+1];
                densidades[3] = m_data[i+1][j+1][k+1];
                densidades[4] = m_data[i]  [j]  [k];
                densidades[5] = m_data[i+1][j]  [k];
                densidades[6] = m_data[i]  [j+1][k];
                densidades[7] = m_data[i+1][j+1][k];

                processaCubo(pontos, densidades);
            }
        }
    }
}

void Grid::processaCubo(QVector3D* pontos, float *densidades) const
{
    QVector3D pontosT[4];
    float densidadesT[4];

    //Triangulacao com cinco tetraedros proposta no livro Fundamentos da Computacao grafica p.266 1a ed.

    //T1
    pontosT[0] = pontos[0]; densidadesT[0] = densidades[0];
    pontosT[1] = pontos[4]; densidadesT[1] = densidades[4];
    pontosT[2] = pontos[1]; densidadesT[2] = densidades[1];
    pontosT[3] = pontos[2]; densidadesT[3] = densidades[2];
    processaTetraedro(pontosT, densidadesT);

    //T2
    pontosT[0] = pontos[5]; densidadesT[0] = densidades[5];
    pontosT[1] = pontos[1]; densidadesT[1] = densidades[1];
    pontosT[2] = pontos[4]; densidadesT[2] = densidades[4];
    pontosT[3] = pontos[7]; densidadesT[3] = densidades[7];
    processaTetraedro(pontosT, densidadesT);

    //T3
    pontosT[0] = pontos[3]; densidadesT[0] = densidades[3];
    pontosT[1] = pontos[2]; densidadesT[1] = densidades[2];
    pontosT[2] = pontos[7]; densidadesT[2] = densidades[7];
    pontosT[3] = pontos[1]; densidadesT[3] = densidades[1];
    processaTetraedro(pontosT, densidadesT);

    //T4
    pontosT[0] = pontos[6]; densidadesT[0] = densidades[6];
    pontosT[1] = pontos[7]; densidadesT[1] = densidades[7];
    pontosT[2] = pontos[2]; densidadesT[2] = densidades[2];
    pontosT[3] = pontos[4]; densidadesT[3] = densidades[4];
    processaTetraedro(pontosT, densidadesT);

    //T5
    pontosT[0] = pontos[1]; densidadesT[0] = densidades[1];
    pontosT[1] = pontos[2]; densidadesT[1] = densidades[2];
    pontosT[2] = pontos[4]; densidadesT[2] = densidades[4];
    pontosT[3] = pontos[7]; densidadesT[3] = densidades[7];
    processaTetraedro(pontosT, densidadesT);

}

QVector3D vertexInterp(float iso, float v1, float v2, const QVector3D& p1, const QVector3D& p2)
{
    float t = (iso - v1)/(v2-v1);
    return (1-t)*p1 + t*p2;
}

QVector3D triangleNorm(QVector3D *tri)
{
    QVector3D v1 = tri[0] - tri[1];
    QVector3D v2 = tri[0] - tri[2];

    return QVector3D::normal(v1,v2);
}

void Grid::processaTetraedro(QVector3D* p, float *d) const
{
    QVector3D tri[2][3];
    QVector3D norm;
    int nt = 0;

    int caso = 0;
    if(d[0] > m_iso)
        caso |= 0x01; //0001
    if(d[1] > m_iso)
        caso |= 0x02; //0010
    if(d[2] > m_iso)
        caso |= 0x04; //0100
    if(d[3] > m_iso)
        caso |= 0x08; //1000

    glColor3f(0.0,0.0,1.0);
    switch(caso)
    {
       case 0x00:
       case 0x0F:
          break;
       case 0x0E:
       case 0x01: // vertice 0 isolado
        //glColor3f(1.0, 0.0, 0.0);
        tri[0][0] = vertexInterp(m_iso,d[0],d[1],p[0],p[1]);
        tri[0][1] = vertexInterp(m_iso,d[0],d[2],p[0],p[2]);
        tri[0][2] = vertexInterp(m_iso,d[0],d[3],p[0],p[3]);
        norm = triangleNorm(tri[0]);
        if (QVector3D::dotProduct(norm, p[0]-p[1]) < 0.0)
            norm = -norm;
        ++nt;
        break;
       case 0x0D:
       case 0x02: // vertice 1 isolado
        //glColor3f(1.0, 0.0, 0.0);
        tri[0][0] = vertexInterp(m_iso,d[1],d[0],p[1],p[0]);
        tri[0][1] = vertexInterp(m_iso,d[1],d[2],p[1],p[2]);
        tri[0][2] = vertexInterp(m_iso,d[1],d[3],p[1],p[3]);
        norm = triangleNorm(tri[0]);
        if (QVector3D::dotProduct(norm, p[1]-p[0]) < 0.0)
            norm = -norm;
        ++nt;
        break;
       case 0x0C:
       case 0x03: // vertices 0 e 1 opostos a 2 e 3

        tri[0][0] = vertexInterp(m_iso,d[1],d[2],p[1],p[2]);
        tri[0][1] = vertexInterp(m_iso,d[0],d[3],p[0],p[3]);
        tri[0][2] = vertexInterp(m_iso,d[0],d[2],p[0],p[2]);

        tri[1][0] = tri[0][0];
        tri[1][1] = tri[0][1];
        tri[1][2] = vertexInterp(m_iso,d[1],d[3],p[1],p[3]);
        norm = triangleNorm(tri[0]);
        if (QVector3D::dotProduct(norm, p[1]-p[3]) < 0.0)
            norm = -norm;
        nt = 2 ;

        break;
       case 0x0B:
       case 0x04: // vertice 2 isolado
        //glColor3f(1.0, 0.0, 0.0);
        tri[0][0] = vertexInterp(m_iso,d[2],d[0],p[2],p[0]);
        tri[0][1] = vertexInterp(m_iso,d[2],d[1],p[2],p[1]);
        tri[0][2] = vertexInterp(m_iso,d[2],d[3],p[2],p[3]);
        norm = triangleNorm(tri[0]);
        if (QVector3D::dotProduct(norm, p[2]-p[0]) < 0.0)
            norm = -norm;
        ++nt;
        break;
       case 0x0A:
       case 0x05: // vertices 0 e 2 opostos a 1 e 3

        tri[0][0] = vertexInterp(m_iso,d[2],d[1],p[2],p[1]);
        tri[0][1] = vertexInterp(m_iso,d[0],d[3],p[0],p[3]);
        tri[0][2] = vertexInterp(m_iso,d[0],d[1],p[0],p[1]);

        tri[1][0] = tri[0][0];
        tri[1][1] = tri[0][1];
        tri[1][2] = vertexInterp(m_iso,d[2],d[3],p[2],p[3]);
        norm = triangleNorm(tri[0]);
        if (QVector3D::dotProduct(norm, p[2]-p[1]) < 0.0)
            norm = -norm;

        nt = 2 ;
        break;
       case 0x09:
       case 0x06: // vertices 0 e 3 opostos a 1 e 2

        tri[0][0] = vertexInterp(m_iso,d[3],d[1],p[3],p[1]);
        tri[0][1] = vertexInterp(m_iso,d[0],d[2],p[0],p[2]);
        tri[0][2] = vertexInterp(m_iso,d[0],d[1],p[0],p[1]);

        tri[1][0] = tri[0][0];
        tri[1][1] = tri[0][1];
        tri[1][2] = vertexInterp(m_iso,d[2],d[3],p[2],p[3]);
        norm = triangleNorm(tri[0]);
        if (QVector3D::dotProduct(norm, p[2]-p[3]) < 0.0)
            norm = -norm;
        nt = 2 ;
        break;
       case 0x07:
       case 0x08: // vertice 3 isolado
        //glColor3f(1.0, 0.0, 0.0);
        tri[0][0] = vertexInterp(m_iso,d[3],d[0],p[3],p[0]);
        tri[0][1] = vertexInterp(m_iso,d[3],d[1],p[3],p[1]);
        tri[0][2] = vertexInterp(m_iso,d[3],d[2],p[3],p[2]);
        norm = triangleNorm(tri[0]);
        if (QVector3D::dotProduct(norm, p[3]-p[0]) < 0.0)
            norm = -norm;
         ++nt;
        break;
    }

    if(nt > 0)
    {
        glBegin(GL_TRIANGLES);
        glNormal3f(norm.x(), norm.y(), norm.z());
        for(int i = 0; i < nt; ++i)
        {
            glVertex3f(tri[i][0].x(),tri[i][0].y(), tri[i][0].z());
            glVertex3f(tri[i][1].x(),tri[i][1].y(), tri[i][1].z());
            glVertex3f(tri[i][2].x(),tri[i][2].y(), tri[i][2].z());
        }

        glEnd();
    }

}
